function stats = BondRetStats(parms,gesol,prob_Nxz,Phi_Nxz_1y,stats_in)

    stats = stats_in;

    % unconditional bond risk premia
    mats = 24:12:60;
    bond_risk_prem_uncond = zeros(size(mats));
    idx_1y = find(gesol.real_bonds.maturities==12);
    p_1y = log(gesol.real_bonds.prices(:,:,:,idx_1y));
    for n = 1:numel(mats)
        idx_buy = find(gesol.real_bonds.maturities==mats(n));
        idx_sell = find(gesol.real_bonds.maturities==mats(n)-12);
        p_buy = log(gesol.real_bonds.prices(:,:,:,idx_buy));
        p_sell = log(gesol.real_bonds.prices(:,:,:,idx_sell));
        bond_risk_prem_uncond(n) = sum(prob_Nxz(:).*(p_sell(:) - p_buy(:) + p_1y(:)));
    end
    stats.bond_risk_prem_real_uncond = bond_risk_prem_uncond;
    
    % conditional bond risk premia, conditional on z
    bond_risk_prem_cond = zeros(numel(mats),parms.Nz);
    
    prob_Nx_given_z = zeros(parms.NN,parms.Nx,parms.Nz);
    for iz = 1:parms.Nz
        prob_Nx_given_z(:,:,iz) = prob_Nxz(:,:,iz)/sum(prob_Nxz(:,:,iz),'all');
    end
    
    for n = 1:numel(mats)
        idx_buy = find(gesol.real_bonds.maturities==mats(n));
        idx_sell = find(gesol.real_bonds.maturities==mats(n)-12);
        p_buy = log(gesol.real_bonds.prices(:,:,:,idx_buy));
        p_sell = log(gesol.real_bonds.prices(:,:,:,idx_sell));
        E_p_sell_1y = reshape(Phi_Nxz_1y*p_sell(:),[parms.NN parms.Nx parms.Nz]);
        for iz = 1:parms.Nz
            bond_risk_prem_cond(n,iz) = sum(prob_Nx_given_z(:,:,iz).*(E_p_sell_1y(:,:,iz) - p_buy(:,:,iz) + p_1y(:,:,iz)),'all');
        end
    end
    
    for iz = 1:parms.Nz
        varname = ['bond_risk_prem_real_z',num2str(iz)];
        stats.(varname) = reshape(bond_risk_prem_cond(:,iz),1,[]);
    end
    
    
    % unconditional volatilities
    bond_ret_uncond_vol = zeros(size(mats));
    bond_exret_uncond_vol = zeros(size(mats));
    
    for n = 1:numel(mats)
        idx_buy = find(gesol.real_bonds.maturities==mats(n));
        idx_sell = find(gesol.real_bonds.maturities==mats(n)-12);
        p_buy = log(gesol.real_bonds.prices(:,:,:,idx_buy));
        p_sell = log(gesol.real_bonds.prices(:,:,:,idx_sell));
        E_p_sell_1y = reshape(Phi_Nxz_1y*p_sell(:),[parms.NN parms.Nx parms.Nz]);
        cond_var_p_sell_1y = reshape(Phi_Nxz_1y*(p_sell(:).^2),[parms.NN parms.Nx parms.Nz]) ...
            - E_p_sell_1y.^2;
        cond_exret = E_p_sell_1y - p_buy + p_1y;
        cond_ret = E_p_sell_1y - p_buy;
        bond_ret_uncond_vol(n) = sqrt(sum(prob_Nxz.*cond_var_p_sell_1y,'all') ...
            + sum(prob_Nxz.*(cond_ret.^2),'all') - sum(prob_Nxz.*cond_ret,'all')^2);
        bond_exret_uncond_vol(n) = sqrt(sum(prob_Nxz.*cond_var_p_sell_1y,'all') ...
            + sum(prob_Nxz.*(cond_exret.^2),'all') - sum(prob_Nxz.*cond_exret,'all')^2);
    end
    
    stats.bond_ret_real_uncond_vol = bond_ret_uncond_vol;
    stats.bond_exret_real_uncond_vol = bond_exret_uncond_vol;
    
end